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1.0 


ABSTRACT 


^  If  internal  wave  measurements  are  to  be  used  for  a  hypothe¬ 
sis  test,  then  the  best  description  of  the  measurements  is  a  multi¬ 
variate  probability  density  function.  This  report  discusses  the 
estimation  of  such  a  probability  distribution,  along  with  the  associ¬ 
ated  concepts  of  Gaussianity,  stationarity,  homogeneity,  isotropy, 
and  frequency- wave  number  spectra  for  random  fields.  By  drawing 
upon  these  concepts,  the  report  offers  suggestions  for  thermistor 
or  current  meter  placement,  for  Doppler  sonar  design,  and  for  tests 
of  Gaussianity  and  stationarity. 
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2.0  INTRODUCTION 

Internal  waves  can  be  described  in  many  different  ways. 

The  description  that  is  most  useful  to  a  detection  theorist  is  a 
statistical  one.  A  detector  is  designed  to  test  two  hypotheses, 

Hj  vs.  Hq.  Under  Hj,  data  consists  of  a  signal  (an  unnatural 
disturbance  of  the  internal  wave  environment),  along  with  naturally- 
occurring  internal  waves  and  noise.  Under  Hq,  data  consists  of 
measurements  of  a  natural  internal  wave  environment,  along  with 
noise  that  is  introduced  by  the  measurement  system.  In  summary 


H^:  Signal  +  natural  background  +  noise 
Hq:  Natural  background  +  noise. 


(1) 


In  order  to  distinguish  between  Hj  and  Hq,  one  can  assign 
costs  to  correct  and  incorrect  judgements,  and  try  to  minimise 
the  expected  risk.  Alternatively,  one  can  fix  the  probability  of  a 
false  alarm  (saying  that  Hj  is  true  when,  in  fact,  HQ  is  true)  and 
find  a  test  that  maximizes  the  probability  of  detection  (correctly 
stating  that  H^  is  true).  In  either  case,  the  resulting  hypothesis 
test  is  implemented  by  forming  a  likelihood  ratio  [1] 


p(data  |Hj) 
p(data  |Hq) 


«1 

^  Threshold 


where 


y 


p(data|H^)  =  the  probability  that  the  measured  data 
would  be  observed  if  H.  were  true, 

i 

i  =  1  or  0. 


(2) 


(3) 
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The  doable  inequality  in  (1)  means  the  detector  says  that  is 
true  if  p(data|Hj)  >  yp(datalHQ),  and  the  detector  decides  that 
Hq  is  true  otherwise. 

From  the  viewpoint  of  a  statistical  hypothesis  test,  one  can 
say  that  the  internal  wave  background  is  specified  when  p(data}H^) 
is  known.  "Environmental  characterization"  is  then  equivalent  to 
finding  p(data  |Hq),  for  any  given  sensor  array.  A  complete  descrip¬ 
tion  of  the  environment  implies  that  the  data  is  measured  in  such  a 
way  that  no  extra  information  would  be  obtained  if  additional  sensors 
were  used. 

This  report  is  concerned  with  measurement  methods  and 
signal  processing  ideas  that  can  be  used  to  obtain  a  complete  speci¬ 
fication  of  the  internal  wave  background,  in  the  sense  that  p(data  Ih^) 
is  estimated,  and  no  additional  sensors  will  change  this  estimate 
(except  to  decrease  measurement  errors). 


3.  0  SAMPLING  INTERVALS  AND  S TA TIONARITY 


3. 1  Sampling  Intervals 

Given  sampled  data  descriptions  over  a  given  interval  of 
time  and  space,  the  number  of  data  samples  depends  upon  the 
sampling  rate.  The  minimum  sampling  rate  that  is  required  to 
completely  specify  a  given  function  is  equal  to  twice  the  bandwidth 
of  the  function  [2]. 

The  temporal  bandwidth  of  internal  waves  is  approximately 
3  x  10  Hz  [3],  so  time  samples  from  any  given  sensor  should  be 
separated  by  less  than  160  sec,  i.  e.  , 

_3 

Temporal  sampling  rate  *  2  "  6  x  10  Hz  (4) 

Sampling  interval  £  (2  Bfc)  1  »  160  sec.  (5) 

The  spatial  bandwidth  can  be  determined  from  wavelength, 

X.  In  the  horizontal  plane  [4,  5] 

20  m  <  X  <  1000  m. 

For  vertically  propagating  internal  wave  components,  the  same 
range  of  wavelengths  is  found  [6].  Thus,  for  both  horizontal  and 
vertical  directions, 

,05  m’1  >  k/2*  >  .001m"1  (6) 

Spatial  sampling  rate  a  2  x  .  05  m  1  (7) 

Sampling  interval  £  10  m.  (8) 
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Equation  (8)  implies  that  the  internal  wave  environment 
should  be  sampled  at  10  m  intervals  (or  less),  in  both  the  horizontal 
and  vertical  directions.  At  each  spatial  sampling  point,  time  sam¬ 
ples  should  be  obtained  at  least  every  160  sec. 

3.  2  Stationarity  Assumptions 

For  any  large  volume  of  ocean,  the  foregoing  analysis  seems 
to  necessitate  the  use  of  three  dimensional  arrays  of  sensors  that 
fill  the  whole  volume,  with  sensors  10  m  apart.  Even  with  Doppler 
sonars,  which  can  remotely  measure  velocity  at  a  given  point,  we 
apparently  need  sonar  beams  that  point  in  all  directions,  in  order  to 
achieve  the  required  sampling  density  in  azimuth  and  elevation. 
Furthermore,  it  is  difficult  to  envision  any  simple  way  to  estimate 
p(data  |Hq),  if  the  joint  probability  of  any  two  measurements  depends 
upon  the  times  at  which  the  measurements  are  taken  (rather  than 
the  difference  in  times). 

If  it  can  be  assumed  that  the  data  is  spatially  and  temporally 
stationary,  then  the  problems  of  dimensionality  and  sampling  are 
considerably  reduced. 

3.  2.  1  Stationarity  in  Time 

Stationarity  in  time  implies  that  p(data  |Hq)  does  not  depend 
upon  the  particular  time  interval  over  which  the  data  is  measured. 
This  assumption,  together  with  ergodicity,  allows  an  observer  to 
estimate  p(datalH^)  by  the  usual  histogram  approach,  or  by  use  of 
Parzen  windows  [7,  8]  or  moment  estimates  [9, 10], 

The  dimensionality  of  the  observation  space  is  still  very 
large,  even  with  a  stationarity  assumption.  Suppose  that  we  meas¬ 
ure  the  responses  of  N  different  spatially  distributed  sensors  at 
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times  t,,  t_,  . . . ,  t. There  are  M  possible  time  differences  in 
this  sequence,  ranging  from  tj-tj  to  t^-t^.  For  each  time  differ¬ 
ence,  there  are  N  different  components  of  the  data.  The  dimension¬ 
ality  of  the  observation  space  is  then  NxM.  For  pdf  estimation,* 
the  number  of  sample  values  should  be  much  larger  than  the  dimen¬ 
sionality  of  the  observation  space. 


3.  2.  2  Spatial  Stationarity 


3.  2.  2.  1  Homogeneity  and  Isotropy 


Two  kinds  of  spatial  stationarity  can  be  defined  [ll].  Given 
a  point  x  =  x^f  x^,  in  a  three  dimensional  space,  we  measure  a 
random  function  a(x)  at  this  point.  The  random  functions  a(x)  for 
various  x-values  or  positions  in  the  space  form  samples  of  a 
random  field. 

The  most  general  version  of  wide-sense  stationarity  for 
such  a  set  of  random  functions  is  to  assume  that  a(x)  is  homogene¬ 
ous.  For  a  homogeneous  random  field  [11], 

E[a(x)3  =  0  (9) 

E  {a(x)  a*(x')3  =  R(x' -x).  (10) 

The  covariance  between  the  random  functions  that  are  measured  at 
different  points  in  the  field  then  depends  only  on  the  magnitude  and 
direction  of  the  vector  that  separates  these  points. 


4 


* 


"pdf"  stands  for  "probability  density  function. 
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A  less  general  version  of  spatial  stationarity  is  a  random 
field  that  is  both  homogeneous  and  isotropic.  For  a  homogeneous, 
isotropic  random  field, 

Efa(x)}  =  0 

E{a(x)a*(x')3  =  R(]]x'-x|]),  (11) 

where  !!x'  -  xl!  is  the  magnitude  of  the  vector  that  separates  x'  and  x. 
An  isotropic  random  field  thus  exhibits  no  direction  dependence,  and 
the  covariance  between  random  functions  measured  at  different  loca¬ 
tions  depends  only  on  the  distance  between  these  locations. 

Both  homogeneous  and  homogeneous -isotropic  random  fields 
have  spectral  distribution  functions  and  corresponding  spectral  density 
functions.  If  the  field  is  not  isotropic,  the  spectral  density  function 
depends  upon  a  vector  wave  number  k: 


S(kj,k2,k2,  co)  =  JJJJ  RfAXj,  Ax2,  Ax3>  At) 


-jfkjAXj  +k2Ax2  +k3Ax3  +  wAt) 


d(Axj)  d(Ax2)  d{Ax3)  d(At), 


where  x'  -  x  has  components  (Ax^,  Ax2,  Ax3). 

For  a  homogeneous-isotropic  random  field 


I C 


S(k,w)  =  //  R(Ax,  At)  e'j(kAx  +  C°At)d(Ax)d(At) 


where 


x'  -  x  =  Ax. 
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3.  2.  2.  2  Implications  and  Applicability  of  Homogeneity 
and  Isotropy  _ _ 

If  the  data  were  both  homogeneous  and  isotropic,  there 
would  be  a  significant  reduction  in  the  required  number  of  sensors 
or  Doppler  sonar  beams.  A  single  beam  or  a  line  of  sensors  would 
be  sufficient  to  generate  the  frequency- wave  number  spectrum  in 
(13). 

The  ocean,  however,  is  not  physically  isotropic  in  three 
dimensions.  In  fact,  the  ocean  is  not  even  homogeneous  in  the 
vertical  direction,  since  many  physical  constants  change  with 
pressure  and  salinity  (hot  to  mention  the  effect  of  surface  and 
bottom).  A  homogeneous -isotropic  assumption  in  three  dimensions 
is  thus  blatantly  unrealistic,  unless  it  is  made  with  respect  to  a 
limited  depth  interval.  Homogeneity  within  such  a  limited  interval 
would  correspond  to  the  concept  of  local  stationarity  [12]. 

Although  it  is  physically  unreasonable  to  assume  that  homo¬ 
geneity  exists  in  the  vertical  direction,  it  would  seem  that  both 
homogeneity  and  isotropy  should  apply  in  the  horizontal  plane,  at 
any  given  depth.  This  observation  implies  that 


E  {a(x)a*(x')3  =  R(Ax12,  x3,  x^). 
In  a  horizontal  plane  at  depth  x^,  where 


we  have 


(15) 
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Sx  (k,  w)  *  JJ  R(Ax12,  x3,  Xj,  At) 

3  -  CO 

-j(k  AXj2  +  wAt) 
e 

d(Axl2)  d(At). 


(17) 


Equation  (17)  implies  that  we  expect  to  see  different  frequency- 
wave  number  spectra  in  horizontal  planes  at  different  depths,  because 
of  the  vertical  inhomogeneity  of  the  ocean. 


to 
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From  our  sampling  results,  Sx3(k,  u>)  should  be  measured 
at  depths  that  are  less  than  ten  meters  apart.  If  a  Doppler  sonar 
beam  is  tilted  in  such  a  way  as  to  measure  internal  wave  velocity 
at  a  range  of  depths  that  exceeds  10  m,  then  we  must  be  careful  in 
our  interpretation  of  the  results.  Unless  it  is  assumed  that  local 
stationarity  exists  within  the  depth  range  that  is  sampled,  the  data 
must  be  segmented  into  10  m  depth  intervals,  and  only  the  data 
from  each  interval  is  then  used  to  estimate  SX3(k,  a>).  Given  this 
analysis  procedure,  a  single  beam  of  a  Doppler  sonar  or  a  single 
line  array  is  sufficient  to  estimate  RfAx^t  x^,  x^). 

The  assumption  of  horizontal  isotropy  can  be  checked  by 
measuring  the  internal  wave  field  along  orthogonal  directions  in  a 
given  horizontal  plane.  If 

R(Axj,  0,  x3,  x3)  i  R(0,  Ax2,  x3,  x^)  (18) 

then  the  process  is  not  isotropic.  In  this  case,  we  have 


(19) 


Sx,(kl'k2'“)  =  JfJ*{Lx V  L*Z'  X3*  x3'  Lt) 

3  -OB 

-jOtjAXj  +  l«2^x2  +  &>At) 
e 

d(AXj)  d(Ax^)  d(At). 

For  a  homogeneous,  non-isotropic  wave  field,  measurements 
become  much  more  difficult.  In  general,  it  is  impossible  to 
measure  R(AXj,  Ax2<  Xy  x^)  with  two  orthogonal  sonar  beams  or 
line  arrays,  since  R(0,  Ax2>  Xy  x3)  and  R(AXj,  0,  x3»  x3)  do  not 
completely  describe  the  autocovariance  function  of  the  data.  For 
homogeneous,  non-isotropic  wave  fields,  a  further  simplifying 
assumption  is  needed  for  convenient  data  gathering. 

3.  2.  2.  3  Separability  for  Homogeneous,  Non-Isotropic 
Random  Fields  _ 

The  most  obvious  simplifying  assumption  for  the  analysis 
of  homogeneous,  non-isotropic  fields  is  separability  of  the  autoco¬ 
variance  function.  By  separability,  we  mean  that 

R(AXj,  Ax2,x3,x3,  At)  ■  R^(AXj,  0,  Xy  x 3#  At) 

(2) 

•  R'  (0,  AX2,  Xy  Xy  At).  (20) 


For  example,  if 


R(Ax,,  Ax,)  =  ZT  c  cosl 
1  2'  mn  \ 


m,  n 


/  2ff  m  Axj ) 

A  T  , 

/  \  T  / 

(21) 
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then  R(Ax^,  Ax^)  is  separable  if 


c  =  a  b  . 
mn  m  n 


In  general,  (20)  implies  that 


(22) 


k2*  w) 


3i1’<kr  “>S*’<V  “»• 

3  3 


(23) 


From  (20),  separability  implies  that  the  random  field  can 
be  completely  described  by  measurements  along  two  orthogonal 
directions,  using  two  orthogonal  sonar  beams  or  line  arrays.  This 
measurement  scenario  is  much  simpler  than  full  two  dimensional 
sampling. 

If  we  plan  to  invoke  a  separability  assumption,  then  we 
should  check  this  assumption  by  measuring  data  along  a  third  direc¬ 
tion  or  at  a  point  outside  the  two  original  lines.  Is  the  covariance 
at  this  point,  relative  to  a  point  at  the  origin,  given  by  the  product 
of  two  other  covariance  functions,  as  in  (20)? 


3.  2.2. 4  Applicability  of  the  Separability  Assumption 

Physically,  it  would  seem  that  separability  will  occur  infre¬ 
quently  in  situations  where  there  is  a  preferred  direction  of  wave 
propagation.  Consider  the  observation  of  wind -propagated  surface 
waves  that  are  all  moving  in  the  same  direction.  The  resulting 
field  may  be  homogeneous,  but  it  is  not  isotropic.  If  the  waves  are 
observed  with  one  axis  pointing  into  the  wind  and  the  other  ortho¬ 
gonal  to  the  direction  of  propagation,  then  the  covariance  of  the 
wave  field  probably  can  be  described  by  a  separable  function.  If, 


however,  the  axes  are  pointing  in  any  other  direction,  then  the 
measured  autocovariance  function  will  not  be  separable.  Separa¬ 
bility  is  thus  orientation-dependent. 

In  situations  where  there  is  no  preferred  direction  of  propa¬ 
gation,  separability  may  be  a  reasonable  assumption.  If  a  check 
shows  that  the  data  is  non- separable,  then  a  rotation  of  coordinates 
may  rectify  the  situation. 
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4. 0  GAUSS1ANITY 

4. 1  Implications  of  Gaussianity 

Homogeneous  random  fields  and  frequency- wave  number 
spectra  are  basic  to  second-order  characterizations  of  internal 
wave  measurements.  If  all  the  data  samples  are  individually  and 
jointly  Gaussian,  then  second-order  characterizations  are  sufficient 
to  completely  describe  the  process. 

4.  2  A  Test  for  Multivariate  Joint  Gaussianity 

We  are  given  a  sequence  of  measurements  a(tj),  .... 

a(t^)  from  N  sensors  at  various  positions  in  the  wave  field,  or  from 
different  resolution  cells  of  a  Doppler  sonar.  We  assume  that 
tn  -  t^  j  »  160  sec,  so  the  samples  are  statistically  independent. 

At  time  tj,  for  example, 

aftj)  =  [ajttj)  a2(t1)  ...  a^)]1,  (24) 

where  aj(t^)  is  the  measurement  from  the  first  sensor  at  time  tj. 

4.  2.  1  The  Standard  Chi-Square  Test 

To  determine  whether  [a,(t  )]  ^  .is  individually  Gaussian, 

i  m  m=  x 

we  can  divide  the  observed  range  of  a^-values  into  I  intervals,  such 
that  there  are  at  least  ten  samples  in  each  interval.  The  mean  and 
variance  of  {aj(tm)}j^_j  are  measured,  and  the  corresponding 
Gaussian  distribution  p(a^)  is  constructed.  The  area  of  the  Gaussian 
distribution  that  is  contained  in  each  of  the  I  intervals  is  determined. 
The  resulting  areas  are  denoted  p^,  and  they  are  the  probability  that 
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a  Gaussian-distributed  data  point  will  fall  within  the  i  interval. 

For  M  data  points,  we  expect  Mp.  points  to  fall  within  the  i**1  interval 

if  the  random  variable  a^  is  Gaussian.  In  fact,  data  points  are 

found  to  fall  within  the  ith  interval,  where  we  have  constructed  the 

interval  size  such  that  V.  >  10. 

1 

A  standard  chi-square  test  with  1-1  degrees  of  freedom  can 
now  be  applied  [13].  Let 

I  (V  -  Mp  )2 

y2  ■  2  (25) 


The  hypothesis  that  the  observed  data  is  Gaussian  is  rejected  if 

X2  >  X2.  (26) 

P 

If  the  observed  data  from  the  first  sensor  really  is  Gaussian,  then 

2  2 

the  probability  that  X  >  Xp  is  p%»  which  is  the  level  of  significance 
of  the  test.  From  Table  1  [13, 14],  if  I  =  5  and  the  level  of  signifi¬ 
cance  is  5%,  we  say  that  the  data  is  non-Gaussian  only  if  >  9.488. 


4.  2.  2  A  Modified  Chi-Square  Test  for  Multivariate  Data 

For  multivariate  data,  N  sensors  or  resolution  cells  form 

an  N-dimensional  observation  space,  and  the  set  of  measurements 

a(t  )  defines  a  sample  point  in  the  space.  The  independent  measure- 

ments  [a(t  )]  ,  yield  M  such  sample  points. 

—  m  m=  l 

The  sample  mean  and  covariance  matrix  are 

t  -  h  2  <"> 


O  RINCON 


■  *  *♦ 


Table  1. 

The  -Distribution  [13, 14] 


Let  p  (x)  be  the  p.  d.  f. ,  the  X  -distribution  with  n  degrees  of  freedom. 
n  2  2 

The  p  percent  value  X  of  X  for  n  d.  of  fr.  is  a  value  such  that  the 

P  2  2 

probability  that  an  observed  value  of  X  exceeds  X  1* 

P 


l-  P 


P  ■  loo"  =  P<X2>X^>  *  /  P.Wdx. 
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940  340  4471  6473  6.30  840  10441  1340  14431  17473  19473  93418  34.70 1140 

3471  440  640  640  7437  9.30  1140  14411  16413  18440  3140  3440  36417  33.30 

4.107  4.70  640  7.01  8.30  9.30  1340  16.113  16.30  19413  3340  36.473  37.30  3440 

440  640  6471  7.70  9437  1040  1840  1640  18.151  3140  3340  3640  29441  3640 

640  640  7431  8447  10487  11.70  1440  1740  19411  33487  34.80  9840  3040  8740 

6413  6414  7.30  9413  1140  13.30  1640  18.413  30.40  23447  3640  99.30  32.30  8940 

6.40  740  840  1040  19.30  13481  1640  19411  31.31*  34.70  3740  3040  33.40  40.70 

7413  7.30  940  10.30  13437  14.40  1740  30.841  38.70  36.80  3840  3240  3440 142411 

740  840  10417  11431  18.n*  1640  1840  3140  33.30  27404  30.10  SS.30  36.10  4340 

840  94 0  10431  19.40  14473  1640  1940  92.77*  36.30  38.411  91.418  36.80  3740  4641* 

840  941*  1140  1940  1640  1740  3040  3340  36471  39.31*  33.871  3640  38.30  46.70 

940  1040  1340  14.841  16414  18.181  2140  24.30  37481  30.30  3340  3740  4040  4840 

1040  1140  1340  1440  17.10  1940  3340  26.813  38.40  33.30  86.172  3840  4140  49.70 

1040  1140  1340  16.40  1840  1940  3340  37.80  39.03  39.10  36.41*  4040  4340  61.10 

1140  1940  14411  16.40  1840  3040  3440  38.10  30.80  3440  3740  4140  44414  6340 

13.10  13.40  16.30  1740  1940  31.70  3640  3940  31.70  3640  8840  4340  46.30  6440 

1340  14.10  16431  18414  30.70  33.717  364  0  304  0  33.813  36.741  40.U3  4440  4640  66.40 

1340  14417  1640  1840  3140  33447  3740  31481  3440  97418  4140  46.0*4840  6640 

1440  16474  17.70  19.70  3347*  9440  9840  23.4*1  3640  3940  43437  4640  4940  6840 

1440  1640  1840  9040  33434  9640  9940  3340  3640  4040  43.70 14740  6040  69.70 
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The  corresponding  N-dimensional  Gaussian  pdf  is 


p(a)  =  [<2ff)N'2lc|1/2]  exp[-j  (a1-^T)g'1(a-^l)]. 


1  ,  T  T^-l. 


"■"I 


jriX-  ^ 

r  C 


ft.  k 


To  generalize  the  chi-square  test,  we  divide  the  observation 

space  into  I  volume  elements  or  bins,  such  that  there  are  at  least 

ten  data  samples  in  each  bin.  The  bins  do  not  necessarily  have  to 

be  the  same  size.  Integrating  p(a)  over  the  i^  volume  element 

results  in  a  probability  p.  that  a  data  sample  will  fall  within  the 
th  * 

i  bin,  if  the  data  is  jointly  Gaussian.  For  M  data  points,  we 

expect  Mp.  points  to  fall  within  the  i^1  bin  if  a  is  jointly  Gaussian. 

In  fact,  »/.  data  points  are  actually  observed  to  fall  within  the  ith 

bin,  where  the  bins  have  been  constructed  such  that  V.  >  10. 

x 

From  this  point  on,  the  chi-square  test  is  identical  to  (25- 
26).  Instead  of  using  the  number  of  samples  that  fall  within  an 
interval  on  the  real  line  for  a^,  we  use  the  number  of  samples  that 
fall  within  an  N-dimensional  volume  element  for  the  lxN  vector  a. 

A  multivariate  Gaussian  pdf  replaces  the  univariate  one  as  a  ref¬ 
erence  function  for  comparison. 


4.  3  What  if  the  Data  Vector  is  Not  Gaussian? 

In  the  vicinity  of  a  source  of  energy  that  is  exciting  internal 
waves,  the  natural  I.  W.  environment  may  give  rise  to  a  non-Gaussian 
data  vector  [15-17].  If  the  chi-square  test  indicates  that  this  is  the 
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case,  then  it  becomes  necessary  to  estimate  the  multivariate  distri¬ 
bution  p(datalH^).  Estimates  of  mean  and  covariance  are  insufficient 
to  specify  a  non-Gaussian  p,  d.  f. ,  and  other  methods  must  be  used. 

Although  the  data  may  be  non-Gaussian  due  to  nonlinear  energy 
transfer  phenomena  [15-17],  it  is  still  reasonable  to  assume  that 
the  data  is  stationary  if  the  energy  source  is  stationary.  For  wind- 
surface,  or  current- seamount  excitation,  a  constant  wind  or  current 
can  be  modelled  as  a  statistically  stationary  source  that  drives  a 
nonlinear,  time-invariant  filter.  The  output  of  this  filter  is  the  I.  W. 
data  at  a  particular  position  relative  to  the  source.  If  the  current  or 
wind  is  a  stationary  process,  then  so  is  the  internal  wave  data  at  a 
fixed  location.  Locations  close  to  the  energy  source  may  yield  non- 
Gaussian  time  series,  whereas  measurements  taken  away  from  the 
source  should  become  Gaussian  [15]. 

Assuming  that  the  data  is  stationary,  the  observation  space 
can  be  divided  into  volume  elements  or  bins,  and  estimates  can  be 
made  of  the  relative  frequency  of  occurrence  of  data  samples  within 
each  element.  Having  obtained  such  a  histrogram,  it  is  tempting 
to  eliminate  discontinuities  by  smoothing  or  interpolating  between 
histogram  values.  For  univariate  probability  distributions,  Parzen 
[7]  has  modelled  the  smoothed  histrogram  estimate  in  terms  of  a 
window  function  that  is  convolved  with  a  sample  distribution.  The 
sample  distribution  is  just  a  sequence  of  delta  functions  scattered 
over  the  observation  space,  where  each  delta  function  represents 
a  measured  data  point.  Meisel  [8]  has  discussed  a  multivariate 
version  of  Parzen' s  estimate.  An  example  of  this  method,  using  a 
Gaussian  window  function,  is  given  in  the  following  paragraphs. 


Let  r  denote  the  random  data  vector,  and  let  a  be  the  m 
—  — m 

sample  vector  or  measurement  of  £.  It  is  convenient  to  use  a  sym¬ 
metric,  separable  window  function 


N 

W(r)  =  n  w(r  ) 
—  n 
n=l 


(30) 


th 


where  r  is  the  random  data  from  the  n  sensor,  representing  the 
th  ® 

n  component  of  the  observation  space. 

To  use  such  a  window  function,  the  data  should  be  normalized 
by  making  the  sample  variances  of  all  observation  space  components 
equal.  The  measured  data  is 


M 

l 

m=  1 


f._)  -  «»ral. 


l  •  •  •  » 


M 


m=l 


(31) 


where  a  is  the  nth  component  of  the  mth  data  measurement.  The 
sample  variance  along  the  n  coordinate  is 


=  -1  2 
n  M  . 

m=l 


M  r 


l  M 

» 

mn  M  ,  ,  in 
i=l 


(32) 


th 


To  normalize  the  data,  all  measurements  from  the  n  sensor  are 

divided  by  9  . 

'  n 

For  a  Gaussian  window  function  w(r  )  with  window  width  9, 

n 

an  estimate  of  the  multivariate  pdf  p(.t!Hq)  is 

M  N 

p(r|H0)  1  ‘  ^  "  ' 


'M' 


m=l  n=  1 


-  a  19  ) 
n  mn  n 


(33) 


.2.N/2-1  M 


=  [M(2ira6’)  ]  S 

m 


A  l  N 

-  exp  |  -  2 
=  1  (  n=  1 


[r  -a  /of/Z9Z\. 
1  n  mn  n  | 
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The  choice  of  window  width,  9,  is  affected  by  several  consid¬ 
erations.  A  wide  window  means  that  many  sample  points  are  averaged 
in  order  to  obtain  p(r|Hg)  at  any  given  point.  From  the  viewpoint  of 
robustness,  this  averaging  would  seem  to  be  desirable.  The  averag¬ 
ing  introduced  by  a  wide  window  will  decrease  the  sensitivity  of  the 
pdf  estimate  to  the  removal  of  any  one  data  point,  and  will  thereby 
increase  the  robustness  of  the  estimate. 

This  reasoning  suggests  that  we  can  allow  the  window  size 
to  become  smaller  as  the  number  of  observations  M  is  increased. 

The  motivation  for  making  the  window  smaller  as  the  data  set  is 
made  larger  follows  from  Parzen's  condition  for  an  unbiased  pdf 
estimate  [7],  This  condition  implies  that  estimates  of  the  form 
(33),  where  w(r)  is  Gaussian  with  standard  deviation  ff(M)  are 
asymptotically  unbiased  if 

lim  ff(M)  =  0.  (34) 

M  —  “ 

The  above  condition,  together  with  the  condition  [7] 

lim  M<T(M)  =  »  (35) 

M  -  - 

assures  that  the  variance  of  the  pdf  estimate  approaches  zero  for  a 
continuous  pdf,  as  M  *•  ®.  In  summary,  (34-35)  are  Parzen’s  condi¬ 
tions  for  unbiased  and  consistent  pdf  estimates.  They  state  that 
the  window  size  should  become  small  as  M  becomes  large,  but  not 
so  small  as  to  violate  (35).  For  any  finite  M,  the  Parzen  pdf  esti¬ 
mates  are  biased  [7], 


k 
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4.  4  What  if  the  Data  Vector  is  Gaussian? 

If  internal  waves  away  from  an  energy  source  result  in 
stationary,  Gaussian-distributed  measurements  in  the  absence  of 
an  unusual  disturbance,  then  an  unusual  disturbance  can  be  detected 
because  it  destroys  the  stationarity  of  the  data.  Haubrich  [18]  has 
devised  a  test  for  stationarity  that  is  based  on  the  fact  that  stationary 
processes  have  uncorrelated  spectral  coefficients  [11].  In  order  to 
distinguish  between  correlated  and  uncorrelated  coefficients,  Hau¬ 
brich  has  measured  the  coherence  between  coefficients  and  has 
compared  the  pdf  of  this  measurement  with  the  pdf  that  is  expected 
for  a  Gaussian,  stationary  process. 

There  are  two  basic  flaws  in  Haubrich' s  approach.  First, 
his  paper  presents  no  independent  test  for  joint  Gaussianity  of  his 
observations.  Thus,  if  the  probability  distribution  of  the  coherence 
estimate  differs  from  the  pdf  corresponding  to  a  stationary,  Gaussian 
I.  W.  process,  it  is  unclear  whether  this  difference  is  caused  by  lack 
of  stationarity  or  lack  of  Gaussianity.  The  same  criticism  applies 
to  bispectral  analysis  [28-34],  A  nonzero  bispectrum  may  be 
caused  by  either  nonlinear  interactions  (which  result  in  non- 
Gaussianity  )  or  by  non- stationarity.  This  defect  is.  remedied  by 
using  a  multivariate  chi-square  test  for  joint  Gaussianity,  as  in 
Section  4.  2.  2.  If  the  data  is  Gaussian,  then  nonzero  coherence 
or  bispectra  must  be  caused  by  nonstationarity. 

The  second  flaw  in  Haubrich’ s  paper  is  that  differences  in 
distributions  are  identified  by  using  histograms  of  "outliers,  "  i.  e. , 
points  that  fall  at  the  extremes  of  the  observation  range.  This 
problem  is  easily  corrected  by  again  using  a  chi-square  test. 
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where  the  reference  distribution  corresponds  to  the  distribution  of 
the  coherence  estimate  for  stationary,  Gaussian  data  [18-20]. 


5.  0  INTERNAL  WAVE  MEASUREMENTS  WITH  DOPPLER 
SONAR 


5.  1  Present  Method 

An  ingenious  method  for  remote  sensing  of  internal  waves 
has  been  conceived  and  implemented  by  Pinkel  [21].  Using  ultra¬ 
sonic  echoes  from  small  volume  scatterers,  Pinkel  can  determine 
the  mean  velocity  of  water  in  a  range-Doppler  resolution  cell.  The 
resulting  plots  of  mean  echo  frequency  as  a  function  of  range  and 
time  illustrate  the  velocity  perturbations  that  are  caused  by  internal 
waves,  and  provide  data  for  construction  of  a  frequency- wave  number 
spectrum. 

In  order  to  simplify  the  signal  processor,  the  Doppler  sonar 
presently  uses  the  echo  covariance  function  to  estimate  the  first 
moment  of  the  echo  power  spectrum  [21,22].  The  concept  of  a 
power  spectrum  implies  stationarity,  whereas  we  actually  expect 
the  center  frequency  to  change  with  time  as  the  wave  moves  through 
a  given  resolution  cell.  The  process  is  thus  based  on  a  local  sta¬ 
tionarity  assumption  [12];  it  is  assumed  that  measurements  can 
be  taken  over  a  temporal  and  spatial  interval  that  is  small  enough 
to  avoid  nonstationarity  effects.  The  use  of  a  first  moment  meas¬ 
ure  also  implicitly  assumes  that  water  motion  is  nearly  uniform 
within  the  resolution  cell  of  the  sonar,  or  at  least  that  the  spectral 
distribution  of  the  echo  is  not  multimodal. 

The  autocovariance  function  of  the  echo  e(t)  is 
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within  any  one  range  resolution  cell.  The  echo  e(t)  has  an  analytic 
representation 


e(t)  =  x(t)+jy(t) 


(37) 


where  y(t)  is  the  Hilbert  transform  of  x(t).  For  any  transmitted 
signal,  a  sample  of  x(t)  is  obtained  by  delaying  the  signal  by  2D/c, 
multiplying  the  echo  data  by  the  signal,  and  integrating  over  the 
duration  of  the  signal.  A  sample  of  y(t)  is  obtained  by  multiplying 
echo  data  by  the  Hilbert  transform  of  the  signal  and  integrating. 
For  sinusoidal  signals,  x(t)  and  y(t)  can  be  obtained  by  quadrature 
demod  ulation. 

mean  frequency  is 

CD 

/  f  S  (f)  df 
0  ee 
e» 

f  S  (f)df 
0 


The 


f  = 


dTRee(T)lT=0  Ree(0) 

2ffj  R  (0)  "  2ffj  R  (0)’ 

““  06 


(38) 


From  (37) 


Ree(T)  =  Ae<T)  exP  [J  Zn  *e(T)l 


where  A  (T)  is  an  even  function  of  T  and  6  (t)  is  an  odd  function  of 
e  e 

T,  since  S  (t)  is  real.  It  follows  that 
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(41) 


Both  the  mean  frequency  and  the  rms  bandwidth  of  the  echo  e(t)  can 
be  determined  from 


Ree(0)  =  E{  |e(t)  |  J 


*  Ree(^)  =  E{e(t)e  (t+A)},  (4: 

i.  e.  ,  from  the  echo  cross-correlation  function  at  a  single  delay  A 
(e.  g.  ,  the  spacing  between  successive  transmitted  pulses)  and  the 
average  power  in  the  echo.  It  would  appear  that  (f  is  not  calculated 
in  the  existing  Doppler  sonar  system. 

5.  2  Possible  Problems  and  Improvements 

5.  2.  1  Spatial  Sampling 

The  1*  beamwidth  of  the  sonar  system  at  75  KHz  means 
that  the  beam  diameter  exceeds  10m  at  ranges  greater  than  600  m. 
Recall  that  10m  was  our  minimum  sample  spacing,  obtained  from 
the  sampling  theorem.  A  10  meter -wide  beam  is  not  the  same  as 
a  10  meter  sampling  interval,  since  spatial  averaging  occurs 
across  the  beam.  Nevertheless,  a  beam  can  be  modelled  as  a 
spatial  impulse  response  function  that  is  convolved  with  the  data. 
The  corresponding  low-pass  spatial  filter  has  a  cut-off  frequency 
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of  sb  0.  05  m  for  a  10  m  wide  beam,  and  this  paesband  should  be 
sufficiently  wide  for  the  data  of  interest,  as  in  (6). 

Pink  el  has  reported  a  range  resolution  of  approximately  25 
m  for  his  system.  Samples  that  are  25  m  apart  are  theoretically 
too  far  apart  for  unique  reconstruction  of  the  environment.  To 
assess  the  practical  significance  of  such  a  large  resolution  cell, 
one  can  try  to  see  whether  Doppler  shifts  sometimes  experience 
a  drastic  change  within  a  single  range  bin.  For  example,  a  very 
high  wave  number  internal  wave  with  a  20  m  wavelength  will  exhibit 
a  large  range  of  velocities  within  a  25  meter  interval. 

The  present  implementation  of  the  Doppler  sonar  assumes 
that  no  appreciable  velocity  changes  occur  over  a  30  m  interval, 
on  the  basis  of  a  coherence  argument  [21].  Having  made  this 
assumption,  only  the  first  moment  of  the  echo  power  spectrum  is 
measured.  This  measurement  method  results  in  a  self-fulfilling 
prophecy;  having  assumed  that  there  is  no  Doppler  spreading  or 
multimodal  spectrum,  we  do  not  see  any. 

In  order  to  check  on  the  sufficiency  of  a  25  m  resolution 
cell,  it  is  necessary  to  refine  the  Doppler  processor  or  to  compare 
identical  data  that  is  measured  with  25  m  and  10  m  range  resolu¬ 
tions. 

To  refine  the  Doppler  processor,  one  can  measure  the 
width  of  the  echo  spectrum  by  using  (41).  This  width,  however, 
may  not  be  very  indicative  of  multiple  velocities  from  a  high- 
wavenumber  event,  since  the  relatively  large  velocities  and 
amplitudes  of  low  frequency-wavenumber  events  will  probably 
dominate  the  spectrum.  It  is  desirable  to  look  at  the  whole  Doppler 


spectrum  from  a  number  of  25  m  range  intervals,  and  to  decide 
whether  such  spectra  may  be  multimodal  or  spread  due  to  high 
frequency-wavenumber  events.  It  should  be  possible  to  perform 
this  analysis  from  tape-recorded  echo  data. 

Another  approach  is  to  measure  the  same  data  with  both 
25  m  and  10  m  resolution  cells.  Better  range  resolution  can  be 
obtained  by  using  a  broader  transmitted  bandwidth.  Bandwidth 
can  be  broadened  without  degrading  the  velocity  estimate  if  coding 
or  modulation  is  used,  i.  e. ,  if  the  signal  has  large  time -bandwidth 
product. 

Another  way  to  improve  the  range  resolution  is  simply  to 
use  a  time  gate  that  isolates  a  10  m  echo  interval  for  each  trans¬ 
mitted  pulse.  If  all  the  scatterers  in  the  10  m  interval  are  moving 
uniformly,  then  they  will  experience  a  range  displacement  equal  to 

8  =  2(v  cos  0  )  A  (44) 

between  two  signals  that  are  A  seconds  apart,  where  v  cos  6  is  the 
radial  velocity  component  of  the  scatterers.  For  a  transmitted 
frequency  f^,  this  displacement  results  in  a  phase  change 

i  -  2(v  cos  0  )  A  (2ff  fg/c).  (45) 

A  sequence  of  pulses  thus  results  in  a  sequence  of  phase  shifts, 
which  show  up  as  a  frequency  shift 

f.  =  2(vcos  0)  f./c  (46) 

d  U 

when  the  echoes  are  Fourier  transformed.  The  accuracy  of  the  f^ 
estimate  depends  upon  the  length  of  the  time  integral  for  calculating 


l  ,iC 

h**J 


the  Fourier  transform,  which  is  proportional  to  the  number  of  gated 
echoes  that  are  coherently  processed. 

In  terms  of  spectral  moment  estimation,  the  gated  echoes 
t  i  N 

en(t)  J  n_j  for  N  transmitted  pulses  with  pulse  repetition  interval 
A  are  used  to  compute 


r(0)  =  s  u  (t)  r 

n=l  n 


R(i)  .  S  .a(t)  Vl<«>.  ( 

n- 1 

It  is  not  necessary  to  store  all  the  pulses  at  the  same  time,  since 


R(0)  =  R^O) 


R(A)  = 


V4> 


where 


R  (0) 
n 


R  (A) 
n 


=  R  (0)  +  |  e  (t)  | 
n- 1  n 


R„-l,i)  +  en(t)  Vl<*>- 


Thus,  only  e^t)  and  ®n  j(t)  need  to  be  stored  at  any  given  time. 
Similarly,  the  Fourier  integral  or  discrete  Fourier  transform  can 
be  divided  into  a  superposition  of  integrals  or  sums  over  disjoint 
time  intervals,  and  this  superposition  can  be  calculated  recursively. 

In  summary,  the  spatial  sampling  that  is  presently  used  by 
die  Doppler  sonar  may  be  too  coarse.  Resolution  can  be  improved 
by  limiting  range  to  less  than  600  m  and  by  using  a  narrow  range 
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gate  with  pulse -to- pulse  coherent  processing,  which  can  be  imple¬ 
mented  recursively. 


5.  2.  2  Increasing  the  Accuracy  of  I.  W.  Velocity  Estimates 

5.2.2. 1  Resolution  Limits 

Internal  waves  with  high  frequencies  and  wavenumbers  tend 
to  propagate  very  slowly.  For  a  two-minute  processing  time, 
velocity  resolution  can  theoretically  be  as  fine  as  1/120  cm/sec 
for  a  75  KHz  center  frequency,  i.  e.  , 

2  fQ  A  v/c  =  1/T 

and 


c  150,000  cm/s  1 

Av  *  ■■■"■  »  — 2-J  ■■■■"  ■ —  =  — -  cm/ sec. 

2f  T  150,000  T  T 


(53) 


The  above  limit  applies  for  an  optimum  (matched  filter)  velocity 
estimator. 


According  to  Pinkel,  such  fine  velocity  resolution  could 
be  wasted,  since  random  motion  of  the  scatterers  will  result  in 
a  minimum  echo  spectral  width  of  1/tc,  where  Tc  is  the  scatter¬ 
ing  correlation  time.  It  would  seem,  however,  that  the  effects 
of  random  scatterer  motion,  along  with  uniform  migration  of 
scatterers,  can  be  significantly  reduced  by  spatial  filtering. 


Spatial  filtering  represents  a  post-processing  method  that 
can  be  used  to  improve  high  resolution  Doppler  data.  The  method 
is  based  on  the  observation  that  the  spatial  frequency  spectrum  of 
scatterer  motion  in  the  absence  of  an  internal  wave  is  usually 
different  from  the  I.W.  spectrum.  Uniform  migration  of  scatterers 
will  result  in  a  Doppler  shift  that  is  nearly  constant  from  one  reso¬ 
lution  cell  to  the  next,  and  the  corresponding  spatial  frequency  is 
near  zero.  Random  motion  of  scatterers  will  usually  be  decorre- 
lated  from  one  resolution  cell  to  the  next,  resulting  in  a  uniform 
distribution  of  spatial  frequencies.  A  spatial  filter  that  passes 
only  the  I.  W.  spectrum  of  interest  (e.  g. ,  high  wavenumber  events) 
should  then  be  capable  of  eliminating  much  of  the  scatterer  motion 
that  is  unrelated  to  the  presence  of  an  internal  wave. 


5.  2.  2.  3  Optimum  Signal  Processing 

The  main  reason  for  the  use  of  covariance  derivatives  to 
generate  spectral  moments  is  the  simplicity  of  the  resulting  sonar 
signal  processor  [21].  Optimum  signal  processing  involves  the  use 
of  matched  filters,  which  can  become  quite  complicated  when  large 
time -bandwidth  signals  are  used  [23].  In  the  near  future,  however, 
fast  Fourier  transforms  will  become  available  as  large-scale  inte¬ 
grated  circuits,  and  a  Fourier  transformer  will  be  just  another 
off-the-shelf  logic  component,  like  a  shift  register,  adder,  or 
accumulator.  Already,  integrated  circuits  exist  for  performing 
many  of  the  operations  in  a  pipeline  FFT  [24].  The  availability 
of  near-real-time,  easy-to-implement  Fourier  transformers  will 
greatly  affect  the  feasibility  of  optimum  signal  processing  techniques 


for  I.  W.  Doppler  sonar.  The  connection  between  Fourier  transforms 
and  velocity  sensing  is  discussed  in  [25,  26].  The  main  points  of  this 
discussion  are  summarized  below. 

The  Doppler  shifts  of  many  different  moving  scatterers  at 
various  ranges  can  be  summarized  in  terms  of  a  scattering  function, 
r(t,  f).  The  scattering  function  describes  the  Doppler  shifts  that 
are  introduced  by  the  medium  as  a  function  of  range.  The  scattering 
function  can  be  measured  by  forming  a  time  window  spectrogram  of 
the  echo  data,  e(t). 

The  spectrogram  is  formed  by  multiplying  the  echo  by  a 
delayed  window  function  w(t  -  t)  and  computing  the  Fourier  transform 
of  the  product: 

CD 

s  (T,d)  *  1  f  e(t)  w(t-T)  exp(-j2ffdt)dt|2.  (54) 

ew  J 

-CD 

Let  w(t)  be  a  time-reversed  version  of  v(t),  i.  e.  , 

v(t)  =  w(-t).  (55) 

Then  [25] 

E[Sev(T,d)}  *  r(t,f)  Vsw(t,f)  (56) 

where  u(t)  is  the  transmitted  sonar  signal  and  ***  denotes  two- 
dimensional  convolution.  In  (56),  T(t,  f)  is  the  scattering  function, 
which  we  would  like  to  measure,  and  Suv(t,f)  is  the  spectrogram 
of  the  transmitted  signal.  E  [  •  }  denotes  ensemble  average. 


Eq.  (56)  says  that  the  echo  spectrogram  ia  a  ameared  ver¬ 
sion  of  the  scattering  function,  where  the  smearing  is  caused  by 
convolution  with  the  apectrogram  of  the  input  signal.  To  reduce  the 
smearing  effect,  the  input  signal  u(t)  and  the  window  function  w(t) 
should  be  designed  to  produce  a  thumbtack-like  signal  spectrogram, 

S  (t,f). 
uv 

The  problem  of  choosing  an  appropriate  signal  u(t)  and 
window  function  w(t)  =  v(-t)  reduces  to  a  standard  signal  design 
problem  by  virtue  of  the  fact  that  [25} 

s„v<T^»  1  lWT-*,|2>  (57> 

where  Xuw*(-T»^)  ia  the  cross -ambiguity  function  of  the  signal  and 
window  functions. 


For  a  narrowband  echo,  an  optimum  estimate  of  center 

frequency  that  is  obtained  from  spectrogram  observations  has  an 

rms  error  <T,  such  that 
a 


a,  ~  [T  T  ] 
d  *•  u  w 


T  *  T 
w  u 


(58) 


[251,  where  T  is  the  signal  duration  and  T  is  the  window  duration. 

L  J  u  w 

The  rms  error  is  minimized  when  the  window  and  signal  have  the 

same  duration,  e.  g. ,  when  the  window  is  matched  to  the  signal. 

* 

In  this  case,  w(t)  =  u  (t),  and 

E[S  (T,rf)}  =  T(t,f)  V  IX,  (-t,f)!2,  (59) 

ev  uu 


where 


Ixuu(-T.rf)!2  ■  I  /  u(t)u*(t-T)exp(-j2*r^t)dt|2 


(60) 


and  Sev(T,  4)  is  given  by  (54-55). 

From  (54),  we  see  that  optimum,  matched-filter  processing 
to  estimate  the  scattering  function  of  the  ocean  (and  hence  to  extract 
Doppler  variations  that  are  induced  by  internal  waves)  is  accomplished 
by  multiplying  echo  data  by  an  appropriate  window  function  and  by 
Fourier  transforming  the  resulting  product.  This  process  should 
become  easier  to  implement  as  LSI  FFT  circuits  become  available, 
or  through  the  use  of  array  processors  [27]  or  pipeline  hardware 
FFT  algorithms  [24]. 


(. 


C 


c 


c 
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5.  2.  2.4  Signal  Design 

If  matched  filtering  or  window  multiplication/FFT  processing 
is  employed,  then  it  is  desirable  to  find  signals  that  have  better  than 
10  m  range  resolution  and  sufficient  Doppler  resolution  for  mapping 
the  slowest  internal  waves  (i.  e.  ,  those  with  high  frequencies  and 
wave  numbers).  For  a  7.  5  m  range  resolution,  the  signal  bandwidth 
should  be  [23] 


B  =  2AR  =  1°°  HZ*  (61) 

For  a  0.  5  cm/s  velocity  resolution,  the  signal  duration  for  f^  =  75 
KHz  should  be  [23] 


T 


2f°4v 


20  sec. 


(62) 


The  required  time -bandwidth  product  is  then  2000.  If  a  coded  wave¬ 
form  is  used,  then  each  pulse  in  the  code  should  have  a  duration  of 
.01  sec,  and  2000  such  pulses  are  needed.  Two  thousand  numbers 
must  be  stored  (e.g. ,  in  a  read-only  memory)  in  order  to  generate 
the  signal  and  the  window  function  or  reference  waveform. 
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The  signal  should  be  designed  to  produce  an  ideal,  thumbtack 
like  ambiguity  function  [23].  Alternatively,  one  can  exploit  the  fact 
that  the  expected  Doppler  shifts  are  limited,  and  the  signal  can  be 
designed  such  that  IXqu(t,  ^)|  =  0  for  [1/20  Hz]  <  |d|  <  the  maxi¬ 

mum  expected  Doppler  shift,  ^max-  A  region  of  the  ambiguity 
plane  would  thus  be  free  of  ambiguity  volume,  producing  a  "volume- 
free  area"  [23].  For  example,  one  could  transmit  a  sum  of  twenty- 
second-long  sinusoids  with  random  center  frequencies,  provided 
that  the  minimum  spacing  between  frequencies  is  greater  than  the 
maximum  possible  Doppler  shift,  as  in  Figure  1.  The  autocorrela¬ 
tion  function  of  the  signal  is  the  inverse  Fourier  transform  of  the 
power  spectrum.  The  use  of  irregularly  spaced  center  frequencies 
should  decrease  the  sidelobe  level  of  the  autocorrelation  function, 
which  would  be  quite  large  for  uniformly  spaced  sinusoids. 


Figure  1.  Power  spectrum  of  a  signal  designed  for  volume-free- 

area  for  Doppler  shifts  between  1/20  Hz  and  d 

rr  max 


5.  2.  3  Identification  of  Sidelobe  Returns 

Pinkel  has  pointed  out  that  a  major  source  of  error  may  be 
echoes  from  scatterers  that  are  not  in  the  main  beam,  i.  e. ,  echoes 
that  enter  the  receiver  through  sidelobes  of  the  beam.  Doppler 
measurements  at  two  different  frequencies  may  help  to  resolve 
this  difficulty. 

If  an  array  of  N  transducer  elements  in  the  x-y  plane  is 
excited  by  a  sinusoid  at  frequency  f,  the  received  signal  at  a  point 
in  the  far  field  that  is  6^  radians  off  boresight  in  the  x  direction  is 

N 

A(0  ,  f)  =  2  exp  [-j2fff(R  -  x  sin0  )/c]  (63) 

n=l  x 

where  R  is  the  distance  from  the  observation  point  to  the  center 

of  the  array,  and  is  the  x- component  of  the  position  of  the  nth 

element,  relative  to  the  array  center.  If  0^  is  less  than  ir/6,  then 

sin  9  m  ft  and 
x  x 

2  N  2 

I A ( f )  ]  =  |  S  exp  [j(2lffxn8inex)/c]| 
n=  1 

N  2 

»  |  w  exp  [j2ffxn(f  Px)/c]  I 
n= 1 

■  lA(f0x)l2.  (64) 

.The  radiation  pattern  of  the  array  is  thus  scaled  when  f  changes. 

The  effect  of  the  scaling  becomes  larger  for  large  angles  away 
from  zero  (boresight).  Near  boresight,  a  frequency  change  has 


35 


ORINCON 


little  effect.  At  points  that  are  removed  from  boresight,  the  radi¬ 
ation  pattern  can  change  considerably. 

We  would  like  to  know  the  contribution  of  scatterers  that 
are  off-boresight  to  the  measured  echo.  This  contribution  can  be 
assessed  by  comparing  I.  W.  echoes  that  are  measured  at  nearly 
the  same  time,  but  at  two  different  frequencies.  A  basic  test  is  to 
cross  correlate  frequency  translated,  "raw"  echo  data.  The  echo 
at  a  given  frequency  is  translated  to  baseband,  using  the  transmitted 
frequency  as  a  local  oscillator  signal.  The  echo  at  a  different  fre¬ 
quency,  measured  at  nearly  the  same  time,  is  processed  in  the 
same  way.  The  two  baseband  echoes  are  random  processes  that 
are  induced  by  scattering  of  a  sinusoid  from  many  randomly- 
distributed  point  targets.  If  the  same  "cloud"  of  point  targets 
gives  rise  to  both  echoes,  the  two  baseband  returns  will  be  highly 
correlated.  If,  however,  different  random  point  targets  are  illumi¬ 
nated,  then  there  should  be  very  little  correlation  between  the  two 
baseband  echoes. 

Sidelobes  that  are  away  from  boresight  will  move  to  new 
angles  when  frequency  is  changed,  and  they  will  illuminate  different 
scatterers  than  before.  The  resulting  baseband  echoes  will  be 
decorrelated  at  the  two  frequencies.  For  signal  energy  at  boresight, 
there  will  be  little  change  in  the  illuminated  scatterers  when  fre¬ 
quency  changes,  and  baseband  echoes  from  scatterers  near  bore¬ 
sight  will  be  highly  correlated  at  the  two  transmitted  frequencies. 

The  correlation  coefficient  between  baseband  echoes  obtained 
at  different  transmitted  frequencies  is  then  a  measure  of  sidelobe 
contributions.  A  correlation  coefficient  that  is  near  unity  indicates 
relatively  small  contributions  from  sidelobes. 


To  identify  Doppler  echoes  from  sidelobes,  one  can  again 
use  the  baseband  echoes  that  are  derived  from  two  different  trans¬ 
mitted  frequencies.  The  power  spectra  of  the  two  echoes  should  be 
different  for  echo  components  that  are  received  through  sidelobes. 
This  difference  should  occur  because  of  a  re -orientation  of  the 
sidelobe  positions  when  frequency  is  changed,  in  accordance  with 
(64).  A  different  angle  between  a  sidelobe  and  the  velocity  vector 
results  in  a  different  Doppler  shift.  Observed  spectral  differences 
in  baseband  echoes  can  be  used  to  construct  a  weighting  function  for 
spectral  conditioning.  Those  Doppler  frequencies  that  exhibit  little 
amplitude  change  are  given  larger  weights  than  Doppler  shifts  that 
exhibit  large  spectral  amplitude  change  when  the  transmitted  fre¬ 
quency  is  changed.  The  resulting  spectral  weighting  function  may 
turn  out  to  be  asymmetric,  which  would  indicate  that  unweighted 
estimates  of  average  Doppler  shift  are  biased. 


5.2.4  Multiple  Beams 

In  Section  3,  the  concepts  of  homogeneity,  isotropy,  and 
separability  were  discussed.  Homogeneity  implies  that  spatial 
autocovariance  functions  depend  only  upon  spatial  shifts  and  not 
upon  the  location  at  which  measurements  are  obtained.  The  shift 
dependence,  however,  can  vary  with  direction  of  shift.  Isotropy 
implies  that  shift  dependence  is  not  direction  dependent;  the  same 
spatial  autocovariance  function  is  measured  no  matter  what  the  shift 
direction.  Separability  means  that  the  autocovariance  function 
R(Ax,  Ay,  Az)  can  be  written  as  a  product  R^(Ax)  R^(Ay)  R^fAz). 

Physically,  homogeneity  is  expected  when  wave  propagation 
occurs  in  a  homogeneous  medium,  without  physical  boundaries  or 


other  spatial  variations.  In  the  ocean,  this  condition  does  not  hold 
in  the  vertical  direction,  but  should  hold  in  the  horizontal  plane  at 
any  given  depth.  Autocovariance  functions  and  frequency- wavenumber 
spectra  must  then  be  indexed  by  a  depth  variable: 


R(Ax,  Ay,  z,  At)  ■  E  {a(x,  y,  z,  t)  A  (x+Ax,  y+Ay,  z,  t+At)J 


and 


Ay,  z.  At) 


-j(k  Ax+k  Ay  +  o>At) 


d(Ax)  d(Ay)  d(At) 


=  S  (k  ,  k  ,  a‘). 
z  x  y 


Isotropy  is  expected  when  there  is  no  preferred  direction  of 
internal  wave  propagation.  This  condition  should  occur  in  a  hori¬ 
zontal  plane  in  deep  water,  away  from  continental  disturbances,  and 
away  from  energy  sources  that  create  internal  waves,  such  as  sea¬ 
mounts.  If  a  random  field  is  isotropic  as  well  as  homogeneous,  then 
it  is  only  necessary  to  measure  the  autocovariance  function  shifts 
in  one  direction,  since  shifts  in  all  other  directions  will  give  the 
same  result.  Thus,  R(Ax,  z.  At)  completely  describes  autocovari¬ 
ance  behavior,  and 

-j(k  Ax+uAt) 

JJJ R(Ax,z,At)e  d(Ax)  d(At) 

=  S  (k  ,  w) 
z  x 


is  a  complete  description  of  the  spectrum  in  a  horizontal  plane  at 
depth  z. 


Although  isotropy  can  often  be  assumed  for  characteriza¬ 
tion  of  the  natural  I.  W.  environment,  it  is  not  a  good  assumption 


if  the  Mine  system  is  to  be  used  for  detecting  unnatural  distur¬ 
bances,  such  as  wakes. 

Separability  is  expected  when  there  is  reason  to  believe 
that  the  autocovariance  function  will  be  symmetric  about  the  x  and 
y  axes.  Symmetry  follows  from  the  fact  that,  if 

R(Ax,  Ay,  z.  At)  =  R  (Ax,  0,  z.  At)  R  (0,  Ay,  z.  At) 
x  y 

then  ]Rx(Ax,  0,  z,  A t)  I  should  be  a  symmetric  function  of  Ax,  in 
order  that  the  power  spectrum  be  real.  Similarly,  |R^(0,  Ay,  z.  At)  | 
is  expected  to  be  a  symmetric  function  of  Ay,  for  any  given  value  of 
Ax. 

On  the  basis  of  the  above  concepts,  some  observations  can 
be  made  about  the  suitability  of  various  beam  patterns  for  Doppler 
sonar  IW  investigations: 

1.  A  vertical  beam  can  be  used  to  assess  the  degree  of  spatial 
inhomogeneity  that  exists  for  internal  waves. 

2.  An  isotropic  assumption  would  allow  data  in  a  horizontal 
plane  to  be  gathered  with  a  single  beam.  At  any  given 
time,  observed  internal  wave  packets  will  probably  be 
non-isotropic.  Isotropy  is  expected  only  for  large  ensemble 
averages  of  IW  wave  measurements.  It  takes  a  whole  day 
for  a  wave  that  moves  at  1  cm/s  to  travel  about  1  km.  To 
gather  enough  data  for  isotropy  to  hold,  one  must  either 
wait  for  many  days,  or  transport  the  sonar  to  a  new  location. 

3.  '  If  two  orthogonal  horizontal  beams  are  used,  then  a 

separability  assumption  is  being  invoked  in  order  to 


reconstruct  a  frequency- wavenumber  spectrum.  Again, 
it  is  doubtful  that  such  an  assumption  holds  for  individual 
wave  packets,  and  it  is  necessary  to  ensemble  average 
measurements  that  are  taken  at  different  locations  or  on 
different  days. 

If  tilted  horizontal  beams  are  used,  then  the  data  from 
different  depths  must  be  processed  separately,  because  of 
vertical  inhomogeneity.  A  beam  that  is  depressed  45*  down 
for  example,  investigates  only  a  30  m  variation  in  Ax,  Ay 
for  each  30-m  thick  horizontal  layer.  Presumably, 

R(Ax,  Ay,  z.  At)  should  be  calculated  for  Ax  and  Ay  greater 
than  30  meters.  A  more  shallow  depression  angle  is  then 
needed,  unless  vertical  homogeneity  holds  over  more  than 
30  m.  An  alternative  is  to  use  many  different  depression 
angles,  which  would  allow  the  sonar  to  gather  data  from 
different  locations  in  a  horizontal  plane  at  each  depth.  Such 
vertical  scanning,  however,  may  be  difficult  to  implement 
with  the  present  system. 


6.  0  SUMMARY  AND  CONCLUSION 


Internal  wave  characterization  has  been  conceived  as  the 
specification  of  a  multivariate  probability  density  function.  The 
number  of  variables  in  this  function  is  determined  by  the  sampling 
theorem,  which  also  indicates  maximum  sensor  spacing. 

Stationarity  assumptions  in  time  and  space  have  been  dis¬ 
cussed,  and  spatial  stationarity  has  been  considered  in  the  form 
of  homogeneity  and  isotropy. 

Autocovariance  functions  and  frequency- wavenumber  spectra 
are  especially  meaningful  if  the  data  is  Gaussian,  as  well  as  stationary. 
A  test  for  multivariate  joint  Gaussianity  has  been  introduced.  If  the 
data  is  not  sufficiently  Gaussian,  one  can  resort  to  estimation  of  the 
probability  distribution  using  a  Parzen  window. 

Doppler  sonar  for  IW  measurements  has  been  discussed,  and 
possible  improvements  have  been  presented.  An  important  aspect 
of  this  discussion  is  the  concept  that  a  Doppler  sonar  should  ideally 
try  to  measure  a  scattering  function,  and  that  this  measurement  is 
best  accomplished  by  using  Fourier  analysis  rather  than  spectral 
moment  estimation.  Even  with  the  present  system,  post-processing 
data  with  a  spatial  filtering  operation  should  help  to  eliminate  the 
effects  of  noise,  random  scatterer  motion,  and  scatterer  migration. 

If  a  Fourier  analysis  capability  is  built  into  the  system, 
then  ROM-generated  window  functions  can  be  used  to  implement 
an  ideal  sonar  processor,  with  better  resolution  and  SNR  than  can 
presently  be  obtained. 


Finally,  the  use  of  multiple  beams  in  a  Doppler  sonar  has 
been  discussed  in  the  context  of  stationarity,  homogeneity,  and 
separability  of  the  random  field  of  internal  wave  measurements. 
Homogeneity  and  separability  assumptions  are  likely  to  hold  only 
for  large  ensemble  averages.  Such  averages  may  suppress  non¬ 


stationary  events  such  as  man-made  wakes.  Estimation  of 
ptdatalH^)  is  thus  much  different  from  estimation  of  p(data|H^). 
It  is  also  likely  that  the  observed  background  for  any  given  man¬ 
made  event  will  not  be  accurately  described  by  p(data|H^)  unless 
the  data  is  observed  over  several  days. 
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